RNA-Seq: differential expression analysis of controls

Libraries required
Load salmon object
Filtering gene counts
- Function
- Counts from all data
Differeitial analysis using limma
vst Counts data and pData
Results
- Save RData files
PCA plots
SessionInfo

Libraries required

library(limma)
library(edgeR)
library(plgINS)
library(sva)
library(SummarizedExperiment)

Load salmon object

load("input/SC_controls_rnaseq_salmon.tds.RData")

Filtering gene counts

Function

counts_filter <- function(counts, nReads = 15, nSamplesPercent = 0.4) {
  counts2 <- counts[apply(counts, 1, FUN = function(x) {
    sum(x >= 15) >= floor(length(x) * 0.4)
  }), ]
  return(counts2)
}

Counts from all data

counts.filt <- counts_filter(counts = salmon@gene.counts)

Differeitial analysis using `limma`

SVA

se <- SummarizedExperiment(assays = counts.filt)
se@colData <- DataFrame(salmon@phenoData)
sv.pl <- svacor(SE = se, form = ~ 0 + Group)

## Using variance-stabilizing transformation

## converting counts to integer mode

## Warning in DESeqDataSet(se, design = design, ignoreRank): some variables in
## design formula are characters, converting to factors

## Number of significant surrogate variables is:  3 
## Iteration (out of 5 ):1  2  3  4  5

limma fit

design <- model.matrix(~ 0 + SV1 + SV2 + SV3 + Group, data = sv.pl@colData)
colnames(design) <- gsub(pattern = "Group", replacement = "", x = colnames(design))

dds <- calcNormFactors(DGEList(counts = counts.filt))
v <- voom(dds, design = design)

contrast.matrix <- makeContrasts(PND15 - PND8, Adult - PND15, levels = design)
fit <- lmFit(v)

fit2 <- contrasts.fit(fit, contrast.matrix)
fit2 <- eBayes(fit2)

DEA results

pnd8.pnd15 <- as.data.frame(topTable(fit2, coef = 1, number = Inf))
pnd8.pnd15 <- data.frame(Genes = rownames(pnd8.pnd15), pnd8.pnd15, stringsAsFactors = F)

pnd15.adult <- as.data.frame(topTable(fit2, coef = 2, number = Inf))
pnd15.adult <- data.frame(Genes = rownames(pnd15.adult), pnd15.adult, stringsAsFactors = F)

Venn Diagram

results <- decideTests(fit2)
vennDiagram(results)

`vst` Counts data and pData

dge <- DGEList(counts = counts.filt, group = salmon@phenoData$Group)
dge <- calcNormFactors(object = dge, method = "TMM")
cpm <- cpm(dge, log = T)
cpm <- removeBatchEffect(x = cpm, covariates = design[, 1:3], design = design[, 4:6])

data <- list(cpm = cpm, vstSV = sv.pl@assays[[2]], voomE = v$E, pData = salmon@phenoData)

save(data,
  file = "./output/data_pData.RData", compress = T,
  compression_level = 3
)

Results

dea.list <- list(
  `PND8 vs PND15` = as.DEA(pnd8.pnd15),
  `PND15 vs Adult` = as.DEA(pnd15.adult)
)

normCounts <- v$E

voomEList <- v

dea.limma <- list(
  `PND8 vs PND15` = pnd8.pnd15,
  `PND15 vs Adult` = pnd15.adult
)

Save RData files

save(dea.list,
  file = "./output/dea_SC_Controls.DEA.RData", compress = T,
  compression_level = 3
)

save(normCounts,
  file = "./output/normCounts_voom_SC_Controls.RData", compress = T,
  compression_level = 3
)

save(voomEList,
  file = "./output/voom_EList_SC_Controls.RData", compress = T,
  compression_level = 3
)

save(dea.limma,
  file = "./output/limma_SC_Controls.RData", compress = T,
  compression_level = 3
)

PCA plots

Raw counts filtered

plPCA(x = counts.filt, samples_data = salmon@phenoData, colorBy = "Group", add.labels = FALSE)

SVA (VST)

plPCA(x = sv.pl@assays[[2]], samples_data = salmon@phenoData, colorBy = "Group", add.labels = FALSE)

CPM

plPCA(x = cpm, samples_data = salmon@phenoData, colorBy = "Group", add.labels = FALSE)

Voom

plPCA(x = v$E, samples_data = salmon@phenoData, colorBy = "Group", add.labels = FALSE)

SessionInfo

devtools::session_info()

## ─ Session info ──────────────────────────────────────────────────────────
##  setting  value                       
##  version  R version 3.6.1 (2019-07-05)
##  os       Ubuntu 16.04.6 LTS          
##  system   x86_64, linux-gnu           
##  ui       X11                         
##  language (EN)                        
##  collate  en_US.UTF-8                 
##  ctype    en_US.UTF-8                 
##  tz       Europe/Zurich               
##  date     2019-09-02                  
## 
## ─ Packages ──────────────────────────────────────────────────────────────
##  package              * version   date       lib source        
##  acepack                1.4.1     2016-10-29 [1] CRAN (R 3.6.1)
##  annotate               1.62.0    2019-05-02 [1] Bioconductor  
##  AnnotationDbi          1.46.1    2019-08-20 [1] Bioconductor  
##  assertthat             0.2.1     2019-03-21 [1] CRAN (R 3.6.1)
##  backports              1.1.4     2019-04-10 [1] CRAN (R 3.6.1)
##  base64enc              0.1-3     2015-07-28 [1] CRAN (R 3.6.1)
##  Biobase              * 2.44.0    2019-05-02 [1] Bioconductor  
##  BiocGenerics         * 0.30.0    2019-05-02 [1] Bioconductor  
##  BiocParallel         * 1.18.1    2019-08-06 [1] Bioconductor  
##  Biostrings             2.52.0    2019-05-02 [1] Bioconductor  
##  bit                    1.1-14    2018-05-29 [1] CRAN (R 3.6.1)
##  bit64                  0.9-7     2017-05-08 [1] CRAN (R 3.6.1)
##  bitops                 1.0-6     2013-08-17 [1] CRAN (R 3.6.1)
##  blob                   1.2.0     2019-07-09 [1] CRAN (R 3.6.1)
##  bookdown               0.13      2019-08-21 [1] CRAN (R 3.6.1)
##  callr                  3.3.1     2019-07-18 [1] CRAN (R 3.6.1)
##  checkmate              1.9.4     2019-07-04 [1] CRAN (R 3.6.1)
##  cli                    1.1.0     2019-03-19 [1] CRAN (R 3.6.1)
##  cluster                2.1.0     2019-06-19 [1] CRAN (R 3.6.1)
##  colorspace             1.4-1     2019-03-18 [1] CRAN (R 3.6.1)
##  crayon                 1.3.4     2017-09-16 [1] CRAN (R 3.6.1)
##  crosstalk              1.0.0     2016-12-21 [1] CRAN (R 3.6.1)
##  data.table             1.12.2    2019-04-07 [1] CRAN (R 3.6.1)
##  DBI                    1.0.0     2018-05-02 [1] CRAN (R 3.6.1)
##  DelayedArray         * 0.10.0    2019-05-02 [1] Bioconductor  
##  desc                   1.2.0     2018-05-01 [1] CRAN (R 3.6.1)
##  DESeq2               * 1.24.0    2019-05-02 [1] Bioconductor  
##  devtools               2.1.0     2019-07-06 [1] CRAN (R 3.6.1)
##  digest                 0.6.20    2019-07-04 [1] CRAN (R 3.6.1)
##  dplyr                  0.8.3     2019-07-04 [1] CRAN (R 3.6.1)
##  edgeR                * 3.26.7    2019-08-13 [1] Bioconductor  
##  evaluate               0.14      2019-05-28 [1] CRAN (R 3.6.1)
##  foreign                0.8-72    2019-08-02 [1] CRAN (R 3.6.1)
##  Formula                1.2-3     2018-05-03 [1] CRAN (R 3.6.1)
##  fs                     1.3.1     2019-05-06 [1] CRAN (R 3.6.1)
##  genefilter           * 1.66.0    2019-05-02 [1] Bioconductor  
##  geneplotter            1.62.0    2019-05-02 [1] Bioconductor  
##  GenomeInfoDb         * 1.20.0    2019-05-02 [1] Bioconductor  
##  GenomeInfoDbData       1.2.1     2019-07-24 [1] Bioconductor  
##  GenomicRanges        * 1.36.0    2019-05-02 [1] Bioconductor  
##  GEOquery               2.52.0    2019-05-02 [1] Bioconductor  
##  ggplot2              * 3.2.1     2019-08-10 [1] CRAN (R 3.6.1)
##  glue                   1.3.1     2019-03-12 [1] CRAN (R 3.6.1)
##  gridExtra              2.3       2017-09-09 [1] CRAN (R 3.6.1)
##  gtable                 0.3.0     2019-03-25 [1] CRAN (R 3.6.1)
##  highr                  0.8       2019-03-20 [1] CRAN (R 3.6.1)
##  Hmisc                  4.2-0     2019-01-26 [1] CRAN (R 3.6.1)
##  hms                    0.5.1     2019-08-23 [1] CRAN (R 3.6.1)
##  htmlTable              1.13.1    2019-01-07 [1] CRAN (R 3.6.1)
##  htmltools              0.3.6     2017-04-28 [1] CRAN (R 3.6.1)
##  htmlwidgets            1.3       2018-09-30 [1] CRAN (R 3.6.1)
##  httpuv                 1.5.1     2019-04-05 [1] CRAN (R 3.6.1)
##  httr                   1.4.1     2019-08-05 [1] CRAN (R 3.6.1)
##  IRanges              * 2.18.2    2019-08-24 [1] Bioconductor  
##  jsonlite               1.6       2018-12-07 [1] CRAN (R 3.6.1)
##  knitr                  1.24      2019-08-08 [1] CRAN (R 3.6.1)
##  later                  0.8.0     2019-02-11 [1] CRAN (R 3.6.1)
##  lattice                0.20-38   2018-11-04 [1] CRAN (R 3.6.1)
##  latticeExtra           0.6-28    2016-02-09 [1] CRAN (R 3.6.1)
##  lazyeval               0.2.2     2019-03-15 [1] CRAN (R 3.6.1)
##  limma                * 3.40.6    2019-07-26 [1] Bioconductor  
##  locfit                 1.5-9.1   2013-04-20 [1] CRAN (R 3.6.1)
##  magrittr               1.5       2014-11-22 [1] CRAN (R 3.6.1)
##  Matrix                 1.2-17    2019-03-22 [1] CRAN (R 3.6.1)
##  matrixStats          * 0.54.0    2018-07-23 [1] CRAN (R 3.6.1)
##  memoise                1.1.0     2017-04-21 [1] CRAN (R 3.6.1)
##  mgcv                 * 1.8-28    2019-03-21 [1] CRAN (R 3.6.1)
##  mime                   0.7       2019-06-11 [1] CRAN (R 3.6.1)
##  miniUI                 0.1.1.1   2018-05-18 [1] CRAN (R 3.6.1)
##  munsell                0.5.0     2018-06-12 [1] CRAN (R 3.6.1)
##  nlme                 * 3.1-141   2019-08-01 [1] CRAN (R 3.6.1)
##  nnet                   7.3-12    2016-02-02 [1] CRAN (R 3.6.1)
##  pillar                 1.4.2     2019-06-29 [1] CRAN (R 3.6.1)
##  pkgbuild               1.0.5     2019-08-26 [1] CRAN (R 3.6.1)
##  pkgconfig              2.0.2     2018-08-16 [1] CRAN (R 3.6.1)
##  pkgload                1.0.2     2018-10-29 [1] CRAN (R 3.6.1)
##  plgINS               * 0.1.5     2019-08-14 [1] local         
##  plotly               * 4.9.0     2019-04-10 [1] CRAN (R 3.6.1)
##  plyr                   1.8.4     2016-06-08 [1] CRAN (R 3.6.1)
##  prettyunits            1.0.2     2015-07-13 [1] CRAN (R 3.6.1)
##  processx               3.4.1     2019-07-18 [1] CRAN (R 3.6.1)
##  promises               1.0.1     2018-04-13 [1] CRAN (R 3.6.1)
##  ps                     1.3.0     2018-12-21 [1] CRAN (R 3.6.1)
##  purrr                  0.3.2     2019-03-15 [1] CRAN (R 3.6.1)
##  questionr              0.7.0     2018-11-26 [1] CRAN (R 3.6.1)
##  R6                     2.4.0     2019-02-14 [1] CRAN (R 3.6.1)
##  RColorBrewer           1.1-2     2014-12-07 [1] CRAN (R 3.6.1)
##  Rcpp                   1.0.2     2019-07-25 [1] CRAN (R 3.6.1)
##  RCurl                  1.95-4.12 2019-03-04 [1] CRAN (R 3.6.1)
##  readr                  1.3.1     2018-12-21 [1] CRAN (R 3.6.1)
##  remotes                2.1.0     2019-06-24 [1] CRAN (R 3.6.1)
##  rlang                  0.4.0     2019-06-25 [1] CRAN (R 3.6.1)
##  rmarkdown              1.15      2019-08-21 [1] CRAN (R 3.6.1)
##  rmdformats             0.3.5     2019-02-19 [1] CRAN (R 3.6.1)
##  rpart                  4.1-15    2019-04-12 [1] CRAN (R 3.6.1)
##  rprojroot              1.3-2     2018-01-03 [1] CRAN (R 3.6.1)
##  RSQLite                2.1.2     2019-07-24 [1] CRAN (R 3.6.1)
##  rstudioapi             0.10      2019-03-19 [1] CRAN (R 3.6.1)
##  S4Vectors            * 0.22.0    2019-05-02 [1] Bioconductor  
##  scales                 1.0.0     2018-08-09 [1] CRAN (R 3.6.1)
##  sessioninfo            1.1.1     2018-11-05 [1] CRAN (R 3.6.1)
##  shiny                  1.3.2     2019-04-22 [1] CRAN (R 3.6.1)
##  SRAdb                  1.46.0    2019-05-02 [1] Bioconductor  
##  stringi                1.4.3     2019-03-12 [1] CRAN (R 3.6.1)
##  stringr                1.4.0     2019-02-10 [1] CRAN (R 3.6.1)
##  SummarizedExperiment * 1.14.1    2019-07-31 [1] Bioconductor  
##  survival               2.44-1.1  2019-04-01 [1] CRAN (R 3.6.1)
##  sva                  * 3.32.1    2019-05-22 [1] Bioconductor  
##  testthat               2.2.1     2019-07-25 [1] CRAN (R 3.6.1)
##  tibble                 2.1.3     2019-06-06 [1] CRAN (R 3.6.1)
##  tidyr                  0.8.3     2019-03-01 [1] CRAN (R 3.6.1)
##  tidyselect             0.2.5     2018-10-11 [1] CRAN (R 3.6.1)
##  usethis                1.5.1     2019-07-04 [1] CRAN (R 3.6.1)
##  vctrs                  0.2.0     2019-07-05 [1] CRAN (R 3.6.1)
##  viridisLite            0.3.0     2018-02-01 [1] CRAN (R 3.6.1)
##  withr                  2.1.2     2018-03-15 [1] CRAN (R 3.6.1)
##  xfun                   0.9       2019-08-21 [1] CRAN (R 3.6.1)
##  XML                    3.98-1.20 2019-06-06 [1] CRAN (R 3.6.1)
##  xml2                   1.2.2     2019-08-09 [1] CRAN (R 3.6.1)
##  xtable                 1.8-4     2019-04-21 [1] CRAN (R 3.6.1)
##  XVector                0.24.0    2019-05-02 [1] Bioconductor  
##  yaml                   2.2.0     2018-07-25 [1] CRAN (R 3.6.1)
##  zeallot                0.1.0     2018-01-28 [1] CRAN (R 3.6.1)
##  zlibbioc               1.30.0    2019-05-02 [1] Bioconductor  
## 
## [1] /home/ubuntu/R/x86_64-pc-linux-gnu-library/3.6
## [2] /usr/local/lib/R/site-library
## [3] /usr/lib/R/site-library
## [4] /usr/lib/R/library

RNA-Seq: differential expression analysis of controls

RNA-Seq: differential expression analysis of controls

Libraries required

Load salmon object

Filtering gene counts

Function

Counts from all data

Differeitial analysis using limma

SVA

limma fit

DEA results

Venn Diagram

vst Counts data and pData

Results

Save RData files

PCA plots

Raw counts filtered

SVA (VST)

CPM

Voom

SessionInfo

Differeitial analysis using `limma`

`vst` Counts data and pData